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Abstract 

We have studied properties of a classical A^-body double-well system coupled to an A#-body 
bath, performing simulations of 2(N$ + Nb) first-order differential equations with N$ — 1 — 10 
and Nb — 1 — 1000. A motion of Brownian particles in the absence of external forces becomes 
chaotic for appropriate model parameters such as Nb, c Q (coupling strength), and {cj n } (oscillator 
frequency of bath): For example, it is chaotic for a small Nb (< 100) but regular for a large Nb 
(> 500). Detailed calculations of the stationary energy distribution of the system fs{u) (u: an 
energy per particle in the system) have shown that its properties are mainly determined by Ng, 
c and T (temperature) but weakly depend on Nb and {co n }. The calculated fs(u) is analyzed 
with the use of the T distribution. Difference and similarity between properties of double-well and 
harmonic-oscillator systems coupled to finite bath are discussed. 
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I. INTRODUCTION 



Many studies have been made with the use of a model describing a classical or quantum 
open system which is coupled to baths consisting of a collection of harmonic oscillators. Such 
a model is conventionally referred to as the Caldeira-Leggett (CL) model pfl E], although 
equivalent models had been proposed earlier by Magalinskii [3] and Ullersma [I]. From 
the CL model, we may derive the Langevin model with dissipation and diffusion (noise) 
terms. Originally the CL model was introduced for iVg-body bath with Nb — > oo, for which 
the Ohmic and Drude-type spectral densities with continuous distributions are adopted. 
Furthermore in the original CL model, the number of particles in a systems, N$, is taken 
to be unity (Ns = 1). We expect that a generic open system may contain any number of 
particles and that a system may be coupled to a bath consisting of finite harmonic oscillators 
in general. In recent years, the CL model has been employed for a study of properties of 
open systems with finite N$ and/or Nb [5Tfl3]. Specific heat anomalies of quantum oscillator 
(system) coupled to finite bath have been studied JSJEIIIS]- A thermalization [HE], energy 
exchange |9J , dissipation [TT] and the Jarzynski equality [131 EH] in classical systems coupled 
to finite bath have been investigated. 

In a previous paper [TU] . we have studied the (Ns + Nb) model for finite iVs-body sys- 
tems coupled to baths consisting of Nb harmonic oscillators. Our study for open harmonic 
oscillator systems with Ns — 1 — 10 and Nb ~ 10 — 1000 has shown that stationary en- 
ergy distribution of the system has a significant and peculiar dependence on Ns, but it 
weakly depends on Nb \1Q\- These studies mentioned above [5TH3] have been made for 
harmonic-oscillator systems with finite Ns and/or Nb- 

Double-well potential models have been employed in a wide range of fields including 
physics, chemistry and biology (for a recent review on double- well system, see Ref. [15]). 
Various phenomena such as the stochastic resonance (SR), tunneling through potential bar- 
rier and thermodynamical properties [H] have been studied. The CL model for the double- 
well systems with Ns = 1 and Nb = oo has been extensively employed for a study on the 
SR [17J. Properties of SR for variations of magnitude of white noise [T7l - [20] and relaxation 
time of colored noise [2TJ E2] have been studied. However, studies for open double-well sys- 
tems with finite Ns and/or Nb have not been reported as far as we are aware of. It would 
be interesting and worthwhile to study open classical double-well systems described by the 
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(Ns + Nb) model with finite Ns and Nb, which is the purpose of the present paper. 

The paper is organized as follows. In Sec. II, we briefly explain the (Ng + Nb) model 
proposed in our previous study [TO]. In Sec. Ill, direct simulations (DSs) of 2(N S + N B ) 
first-order differential equations for the adopted model have been performed. Dynamics of 
a single double- well system (Ng = 1) coupled to a finite bath (2 < Nb < 1000) in the phase 
space is investigated (Sec. Ill B). We study stationary energy distributions in the system 
and bath, performing detailed DS calculations, changing N s , N B , the coupling strength and 
the distribution of bath oscillators (Sec. Ill C). Stationary energy and position distributions 
obtained by DSs are analyzed in Sec. IV. The final Sec. V is devoted to our conclusion. 



II. ADOPTED (N s + N B ) MODEL 



We consider a system including Ng Brownian particles coupled to a bath consisting of 
independent Nb harmonic oscillators. We assume that the total Hamiltonian is given by 



H = H s + H B + Hi, 
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(3) 
(4) 



Hi = - c kn (Qk - q n ) 2 , 

k=l n=l 

where Hs, Hb and Hi express Hamiltonians for the system, bath and interaction, respec- 
tively. Here M (m) denotes the mass, Pk (p n ) the momentum, Qk (q n ) position of the 
oscillator in the system (bath), V(Qk) signifies the potential in the system, u n stands for 
oscillator frequency in the bath, and c n k is coupling constant. The model is symmetric with 
respect to an exchange of system -H- bath if V(Q) is the harmonic potential. From Eqs. 
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Q-Q, we obtain 2(Ns + Nb) first-order differential equations, 



ft = g, (5) 

N B 

Pk = -V'(Qk) -^c fen ((5 fc - q n ), (6) 

n=l 

In = — , (7) 

m 

N s 

p n = -mu 2 n q n - c kn (q n -Q k ), (8) 
fe=i 



which yield 



N B 

MQ k = -V'{Q k ) - cUQk - q n ), (9) 

n=l 

N S 

mq n = -mu 2 n q n - } j c kn (q n - Q k ), (10) 

k=l 

with prime (') and dot (■) denoting derivatives with respect to the argument and time, 
respectively. It is noted that the second term of Eq. ^ or ^ given by 

N B 

F ( k eff) = -J^c kn (Qk-q n ), (11) 

n=l 

plays a role of the effective force to the fcth system. 
A formal solution of Eq. (10) for q n (t) is given by 

q n {t) = gn(0) cos cl> w t + ^ w j^ sin u n t + V] — ^- / sin a) n (t - t')Qe(t') dt', (12) 

with 

W n = ^ 7 = W n + 7 • ( 13 ) 

k=l k=l 

Substituting Eq. (12) to Eq. (|9|, we obtain the non-Markovian Langevin equation given by 

N s N s » t 

MQ k (t) = -V'(Q k ) - Mj2tkeQe(t) -J2 / lu{t-t')Q t {t') dt' 

l=\ i=i ^° 

N s 

- ^lht{t)Q t (0) + (k(t) {k = 1 to N s ), (14) 
e=i 
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with 



Ck(t) 



N B 
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COS UJ n t, 



q n (u) cos w n t H — sinwj 



0J r , 



(15) 
(16) 
(17) 



where denotes the additional interaction between and £th particles in the system 
induced by couplings {ck n }, 1ke(t) the memory kernel and (k the stochastic force. 
If the equipartition relation is realized in initial values of q n (0) and q(0), 

(mu 2 n q n (0) 2 ) B = (mq n (0) 2 ) B = k B T, (18) 
we obtain the fluctuation-dissipation relation: 

(Cfc(*)Ck(0>B = k B T lkk (t-t'), (19) 
where (-) B stands for the average over variables in the bath. 



In the case of N B — > oo, summations in Eqs. ( 15 )-( 17) are replaced by integrals. When 
the spectral density defined by 



J(u) 



V — ^5(w - u„ 



is given by the Ohmic form: J(u) oc u for < w < Wd, the kernel becomes 

sinw D t 
7(t) cx — oc b(t), 

Tit 

which leads to the Markovian Langevin equation. 



(20) 



(21) 



In the case of N$ = 1, we obtain £ and 7 in Eqs. (15) and (16) where the subscripts k 
and £ are dropped (e.g., Ck n = c n ), 

M£(t) =VJl-^ , (22) 



7(0 
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The additional interaction vanishes (£ = 0) if we choose c n = mu^ in Eq. (22). 
In the case of Ns 7^ 1, however, it is impossible to choose {ck n } such that 



(23) 



is 



realized for all pairs of (k, £) in Eq. (15). Then Q k is inevitably coupled to Qi for £ 7^ k with 
the superexchange-type interaction of antiferromagnets: — Y^ n Cfc n Q n /mu)^ in Eq. (15). 



III. MODEL CALCULATIONS FOR DOUBLE- WELL SYSTEMS 



Calculation methods 



2\2 



(24) 



We consider a system with the double-well potential 

V(Q) = (m) iQ2 - Q 

which has the stable minima of V(±Qo) = at Q = ±<5o and locally unstable maximum of 
V(0) = A at Q = with the barrier height A. We have adopted Qo = 1.0 and A = 1.0 in 
our DSs. 

It is easier to solve 2(Ns + Nb) first-order differential equations given by Eqs. Q-Q than 



to solve the N s Langevin equations given by Eqs. (14)- (17) although the latter provides us 



with clearer physical insight than the former. In order to study the Ns and Ng dependences 
of various physical quantities, we have assumed that the coupling Ck n is given by pm El] 

Ckn = iv~^' (25) 

because the interaction term includes summations of J2k=i an< ^ ^n=i m -^q- If is noted 
that with our choice of Ck n , the interaction contribution is finite even in the thermodynamical 
limit of Nb — > oo because the summation over n runs from 1 to Nb in Eq. Q. DSs of Eqs. 
(J5|-(j8]) have been performed with the use of the fourth-order Runge-Kutta method with 
the time step of 0.01. We have adopted ks = 1.0, M = m = 1.0, c Q = 1.0, and uj n = 1.0 
otherwise noticed. 

We consider energies per particle u v (t) in the system (^=S) and the bath (r/=B) which 
are assume to be given by 

N S 
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(26) 
(27) 



which is valid for the weak interaction, although a treatment of the finite interaction is 
ambiguous and controversial [5J E] . 
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B. Dynamics of a particle in the (Q,P) phase space 



1. Effect of c 

First we consider an isolated double-well system {N$ — 1 an d c G = 0.0). Figure [l] shows 
the phase-space trajectories in the (Q, P) phase space for this system with six different initial 
system energies Es - For Es = 0.0, the system has two stable fixed points at (Q, P) = 
(±1.0,0.0), and for E So = A = 1.0 it has one unstable fixed point at (Q,P) = (0.0,0.0). 
In the case of 0.0 < Eg < 1.0, the trajectory is restricted in the region of Q > 0.0 (or 
Q < 0.0). In contrast in the case of Es > 1-0, trajectory may visit both regions of Q > 0.0 
and Q < 0.0. The case of Ego = 1-0 is critical between the two cases. 

Next the double-well system is coupled to a bath. In our DSs, we have assumed that 
system and bath are decoupled at t < where they are in equilibrium states with Es = T, 
the temperature T being defined by T = ub- We have chosen initial values of Q(0) = 1.0 



and P(0) = y/2M[E So - V(Q(0))] for a given initial system energy Eg - Initial conditions 
for q n (0) and p n (0) are given by random Gaussian variables with zero means and variance 



proportional to T [Eq. (18)] [10J. Results to be reported in this subsection have been 
obtained by single runs for t — to 1000. 

Figures |2]( a) and|2](b) show a strobe plot in the (Q, P) phase space (with a time interval 
of 1.0) and the time-dependence of Q(t), respectively, for E So = 1.0, N s = 1, N B = 100 and 
c Q = 0.2. The trajectory starting from Q(0) = 1.0 goes to the negative-Q region because a 
particle may go over the potential barrier with a help of a force (noise) originating from bath 
given by Eq. (11). The system energy fluctuates as shown in Fig. [21(c), whose distribution 



is plotted in Fig. |2| 

Results in Fig. [2] are regular. In contrast, when a coupling strength is increased to 
c a = 1.0, the system becomes chaotic as shown in Figs. [3]^a) and[3^b) where a strobe plot 
in the (Q, P) phase space and the time-dependence of Q(t) are plotted, respectively. This is 
essentially the force- induced chaos in classical double- well system [23]: although an external 
force is not applied to our system, a force arising from a coupling with bath given by Eq. 



( 11 ) plays a role of an effective external force for the system. Figures |3](c) and[3[d) show that 
in the case of c Q = 1.0, u$ has more appreciable temporal fluctuations with a wider energy 
distribution in fs(u) than in the case of c a = 0.2. Although system energies fluctuate, they 
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FIG. 1: Plot of phase-space trajectories for a particle in an isolated double-well system (c = 0.0). 
Trajectories are plotted for energies of Es / 'A = 0.0, 0.5, 0.8, 1.0, 1.2 and 1.5. 

are not dissipative at 0.0 < t < 1000.0 in DSs both for c Q = 0.2 and c Q = 1.0 with N B = 100. 

2. Effect of L0 n distributions 

We have so far assumed to n — 1.0 in the bath, which is now changed. Figures [4ta) and 
|2tc) show strobe plots for co n = 0.5 and 2.0, respectively, which are regular and which are 
different from a chaotic result for uj n = 1.0 shown in Fig. |4^b). When we adopt {co n } which 
is randomly distributed in [0.5, 2.0], a motion of a system particle becomes chaotic as shown 
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FIG. 2: (Color online) (a) Strobe plot in the (Q,P) phase space (with a time interval of 1.0), 
(b) Q(t), (c) us(t), and (d) the system energy distribution fs(u) obtained by a single run for 
E So = 1-0, N s = 1, N B = 100, T = 1.0 and c Q = 0.2. 

in Fig. |4^d). This is because contributions from oj n ~ 1.0 among [0.5.2.0] induce chaotic 
behavior. 
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FIG. 3: (Color online) (a) Strobe plot in the (Q,P) phase space, (b) Q(t), (c) us(t), and (d) 
the system energy distribution fs(u) obtained by a single run for Eg a = 1.0, N$ = 1, Nb = 100, 
T = 1.0 and c Q = 1.0. 

3. Effect of N B 

We have repeated calculations by changing Nb, whose results are plotted in Figs. |5](a)- 
^d). Figures [5^ a) , Mh) and[5^c) show that chaotic behaviors for Nb = 2 and Nb = 10 are 
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FIG. 4: Strobe plots in the (Q, P) phase space for various distribution of {oj n }: (a) u n = 0.5, (b) 
D n = 1.0, (c) u n = 2.0 and (d) w„ E [0.5,2.0] obtained by single runs with Ego = 1-0, N$ = 1, 
iVij = 100, T = 1.0 and c G = 1.0. 

more significant than that for iVg = 100. On the contrary, chaotic behavior is not realized 
for Nb = 1000 in Fig. |5^d), which is consistent with the fact that chaos has not been 
reported for the double-well system subjected to infinite bath. 
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FIG. 5: Strobe plots in the (Q,P) phase space for various Nb'- (a) Nb = 2, (b) Nb = 10, (c) 
N B = 100 and (d) N B = 1000 with E So = 1.0, N S = 1, T = 1.0 and c a = 1.0. 

4- Effect of initial system energy E$ 

Next we change the initial system energy of Es - Figures |6^a)-(d) show strobe plots in 
the (Q,P) phase space for various Es with N$ = 1, Nb = 100, T = 1.0 and c Q = 1.0. 
Figure [6]^a) shows that for Es Q = 0.5, the regular trajectory starting from Q = 1.0 remains 
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in the positive-Q region because a particle cannot go over the potential barrier of A = 1.0. 
For Eso = 0.8, chaotic trajectories may go to the negative-Q region with a help of force from 
bath [Eq. [II]. Figure [6^d) shows that when E So is too large compared to A (E So / A = 1.2), 
the trajectory again becomes regular, going between positive- and negative-Q regions. 

Figure [7] shows the system energy distribution fs(u) for various Ego- fs(u) moves upward 
as Eso is increased. It is noted that peak positions of fs(u) for Es = 0.5 — 1.0 locate at 
u ~ 1.0 while that for E$o = 1-2 locates at u ~ 1.35. 



C. Stationary energy probability distributions 



In this subsection, we will study stationary energy probability distributions of system and 
bath which are averaged over N r (=10 000) runs stating from different initial conditions. 
Assuming that the system and bath are in the equilibrium states with T = ub = u$ at 
t < 0.0, we first generate exponential derivatives of initial system energies {Ej}: p{Ej) oc 
exp(— (3Ej) (J = 1 to NsN r ) for our DSs where (3 = 1/ksT. A pair of initial values of Qj(0) 
and -Pj(O) for a given Ej is randomly chosen such that they meet the condition given by 
Ej = Pj(0) 2 /2M + V(Qj(0)). The procedure for choosing initial values of q n (0) and p n (0) 
is the same as that adopted in the preceding subsection [10]. We have discarded results for 
t < 200 in our DSs performed for t — to 1000. 

Before discussing cases where N$ and Nb may be greater than unity, we first study a 
pedagogical simple case of N$ = Nb = 1: a particle with double-well potential is subjected 
to a single harmonic oscillator. Double-chain curves in Fig. |8](a) and [81(b) show energy 
distributions of the system [f s (u)} and bath [f B (u)}, respectively, with c D = 0.0, where 
u = Ug (u = Ub) for the system (bath). Both fs{us) and follow the exponential 

distribution because the assumed initial equilibrium states of decoupled system and bath 
persist at t > 0.0. When they are coupled by a weak coupling of c Q = 0.1 at t > 0.0, fs(u) 
and /b(m) almost remain exponential distributions except for that fs(u) has a small peak 
at it = 1.0, as shown by dashed curve in Fig. [8]^a). This peak has been realized in Figs. 
|2]^d) and[3^d). It is due to the presence of a potential barrier with A = 1.0 in double- well 
potential because the peak at u = 1.0 in fs{u) is realized even when T 7^ 1.0, as will be 



discussed later in 4- Effect ofT (Fig. 12). This peak is developed for stronger couplings of 
c = 1.0 and 2.0, for which magnitudes of fs(u) at small u are decreased, as shown by solid 
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FIG. 6: Strobe plots in the (Q,P) phase space for various Ego- (a) Es = 0.5, (b) Ego = 0.8 (c) 
E So = 1.0 and E So = 1.2 with iV 5 _ lj jv" b = 100, T = 1.0 and c = 1.0. 

and chain curves in Figs. (8^a) and|8](b). 
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FIG. 7: (Color online) System energy distributions fs(u) for Ego = 0.5, 0.8, 1.0 and 1.2 with 
Ns = 1, Nb = 100, T = 1.0 and c Q = 1.0, curves being successively shifted upward by two for 
clarity of figures. 

1. Effect of c 

We change the coupling strength of c a . Figures [9^ a) and [9]^b) show fs(u) and /b(w), 
respectively, for c D = 0.2, 1.0, 5.0 and 10.0 with iV 5 = 1, N B = 100 and T = 1.0. / 5 (u) for 
c a = 0.2 nearly follows the exponential distribution. When c Q becomes larger, magnitudes 
of fs(u) at u < 1.0 are decreased while that at u > 1.0 is increased. In particular, the 
magnitude of /s(0) is more decreased for larger c Q . 
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FIG. 8: (Color online) Stationary distributions of (a) fs(u) and (b) for various c a : c Q = 0.0 

(double-chain curves), 0.1 (dashed curves), 1.0 (solid curves) and 2.0 (chain curves) obtained by 
10 000 runs with N s = N B = 1 and T = 1.0. 

2. Effect of u n distributions 

Although we have assumed oj n — 1.0 in bath oscillators, we will examine the effect of their 
distribution, taking into account two kinds of random distributions given by u n G [0.5, 2.0] 
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FIG. 9: (Color online) Stationary distributions of (a) fs(u) and (b) for various c a : c Q = 0.2 

(dashed curves), 1.0 (chain curves), 5.0 (double-chain curves) and 10.0 (solid curves) with Ns = 1, 
N B = 100 and T = 1.0. 



and u n G [2.0,3.0]. From calculated results shown in Figs. [I0|(a) and 10 'b), we note that 
fsiu) and /b('u) are not much sensitive to the distribution of {co n } in accordance with our 
previous calculation for harmonic oscillator system [TOl 125] . This conclusion, however, might 
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not be applied to the case of infinite bath where distribution of {u n } becomes continuous 
distribution. Ref. [H] reported that the relative position between oscillating frequency ranges 
of system and bath is very important for a thermalization of the harmonic oscillator system 
subjected to finite bath. 

3. Effect of N B 

We have calculated fs(u) and /b('u), changing Nb but with fixed Ns = 1, whose results 
are shown in Figs. [TT|^a) and[PTjfo). For larger Nb, the width of /b(u) becomes narrower as 
expected. However, shapes of fs(u) are nearly unchanged for all cases of Nb = 1, 10, 100 
and 1000. 



I Effect of T 

We change the temperature of the bath. Figures |l2|a) and[l~2"|(b) show fs(u) and /b(u), 
respectively, for T = 0.5, 1.0 and 1.5 with Ng = 1, Nb = 100 and c Q = 1.0. When T 
is decreased (increased), positions of /b('u) move to lower (higher) energy such that mean 
values of ub correspond to T. For a lower temperature of T = 0.5, magnitude of fs(u) at 
u < 1.0 is increased while that at u > 1.0 is decreased. The reverse is realized for higher 
temperature of T = 1.5. We should note that the peak position in fs{u) at u — 1.0 is not 
changed even if T is changed because this peak is related to the barrier with A = 1.0 of the 
double-well potential. 



5. Effect of N s 

Although N s = 1 has been adopted so far, we will change N s to investigate its effects on 



stationary energy distributions. Figure 13 ^a) shows fs(u) for Ns = 1, 2, 5 and 10. fs(u) 
for Ns = 1 shows an exponential-like distribution with fs{0) ^ at u = 0.0. In contrast, 
fs(u) vanishes at u = 0.0 for Ns = 2, 5 and 10. Figure |T3~)^a) shows that shapes of fs(u) 



much depend on Ns while those of /b(u) are almost unchanged in Fig. [l3|b). 
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FIG. 10: (Color online) Stationary distributions of (a) fs(u) and (b) /b(w) for various distributions 
of {uj n }: uj = 1.0 (dashed curves), uj n S [0.5,2.0] (solid curves) and ui ri G [2.0,3.0] (chain curves) 
with N s = 1, N B = 100, T = 1.0 and c c = 1.0. 
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FIG. 11: (Color online) Stationary distributions of (a) fs(u) and (b) /b(w) for various Nb- 
Nb = 1 (solid curves), 10 (dashed curves), 100 (chain curves) and 1000 (bold solid curves) with 
N s = 1, T = 1.0 and c a = 1.0. f B (u) for Nb = 1000 is multiplied by a factor of 1/3. 
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FIG. 12: (Color online) Stationary distributions of (a) fs(u) and (b) for various T: T = 0.5 

(chain curves), 1.0 (solid curves) and 1.5 (dashed curves) with N$ = 1, Nb = 100 and c Q = 1.0. 
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FIG. 13: (Color online) Stationary distributions of (a) fs(u) and (b) /b(w) for various Ns- N$ = 1 
(dashed curves), 2 (doted curves), 5 (chain curves) and 10 (solid curves) with Nb = 100, T = 1.0 
and c = 1.0. 
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IV. DISCUSSION 



A. Analysis of stationary energy distributions 

Our DSs in the preceding section have shown that fs(u) depends mainly on Ns, c Q and 
T while /b(w) depends mostly on Nb and T for Ns <C Nb- We will try to analyze fs(u) 
and /b(m) in this subsection. It is well known that when variables of Xi (i = 1 — iV) are 
independent and follow the exponential distributions with the same mean, the distribution of 
its sum: X = £\ is given by the T distribution. Then for an uncoupled system (c = 0.0), 
fs{u) and are expressed by the T distribution given by [TU] 

1 

with 



/„(«) = u a ^ l e-^ u = g{u;a n ,b v ), (28) 

Zlr, 



a v = N n , b v = N v p, (29) 
Z v = ( 30 ) 

Or) 

where 77 =S and B for a system and bath, respectively, and T(x) is the gamma function. In 
the limit of Ns = 1, the T distribution reduces to the exponential distribution. Mean (fj, v ) 
and variance (<r^) of the T distribution are given by 

= T~> ^ = W, (31) 

from which and 6,, are expressed in terms of n n and a r] 

a, = §, \ = (32) 

We have tried to evaluate fs(u) and Jb{u) for the coupled system (c G 7^ 0.0) as follows: 
From mean (//,,) and root-mean-square (RMS) (a v ) calculated by DSs, and 6^ are deter- 



mined by Eq. (32), with which we obtain the T distributions for fs(u) and /#(«)• Filled 



and open squares in Fig. 14 show \ib and as, respectively, as a function of N s . We obain 



Hb = 1-0 and <7b = 0.1 nearly independently of Ns, which yield as = &b = 100.0 in agree- 



ment with Eq. (29). Filled and open triangles in Fig. 14 express the Ns dependence of fj,s 
and as obtained by DSs with c a = 1.0, Nb = 100 and T = 1.0. Calculated mean and RMS 
values of (n 8 , cr s ) are (1.07,0.98), (0.99,0.70), (0.99,0.44) and (0.99,0.319) for N s = 1, 2, 5 

23 
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FIG. 14: (Color online) Ns dependences of fi v and a v of system (77 = S) and bath (77 = -B) with 
iVe = 100 and T = 1.0: filled (open) triangles denote fjL$ ( a s) with c D = 1.0: filled (open) circles 
express fj,s (as) with c Q = 10.0: filled (open) squares show /ig (&b) with c Q = 1.0. 



and 10, respectively, for which Eq. (Q yields (a s , 65) = (1.18, 1.11), (2.04, 2.05), (4.50, 5.06) 
and (9.82, 9.97). These values of as and 65 are not so different from Ns and Ns(3 given by 
Eq. ((29l. We have employed the T distribution with these parameters as and bs for our 



analysis of fs(u) having been shown in Fig. 13 a). Dashed curves in Figs. [T5|(a)-(d) express 
calculated V distributions, which are in fairly good agreement with fs(u) plotted by solid 
curves, except for Ns = 1 for which g(0) = 0.0 because as = 1.18 > 1.0 while /s(0) 7^ 0.0. 
Similar analysis has been made for another result obtained with a larger c Q = 10.0 for 
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FIG. 15: (Color online) u dependences of f s (u) for (a) N s = 1, (b) N s = 2, (c) N s = 5 and (d) 
Ns = 10 with T = 1.0, c = 1.0 and iV^ = 100 obtained by DSs (solid curves): dashed and chain 
curves express V and exponential distributions, respectively (see text). 

Nb = 100 and T = 1.0. ^-dependences of calculated fis and as are plotted by filled 



and open circles, respectively, in Fig. 14 Calculated (/j,s,<Js) are (2.88,2.61), (1.81,1.12), 



;i.21,0.47) and (1.07,0.31) for N s = 1, 2, 5 and 10, respectively, which lead to (a s ,b s ) 
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FIG. 16: (Color online) u dependences of f s (u) for (a) N s = 1, (b) N s = 2, (c) N s = 5 and (d) 
Ns = 10 with T = 1.0, c = 10.0 and iV^ = 100 obtained by DSs (solid curves): dashed and chain 
curves express T and exponential distributions, respectively (see text). 



1.22,0.42), (2.62,1.45), (6.57,5.44) and (11.54,10.81) by Eq. (32). Obtained a s and b s are 



rather different from Ns and Ns(3 given by Eq. (29). Dashed curves in Figs. b6[a)-(d) show 



T distributions with these parameters, which may approximately explain fs(u) obtained by 
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DSs in the phenomenologically sense, except for Ns = 1 for which g(0) = 0.0 but /s'(O) 7^ 0.0. 



We note in Fig. 15 a) or [I6[a) that an agreement between g{u) and fs(u) with Ns = 1 
is not satisfactory. We have tried to obtain a better fit between them, by using the q-T 
distribution g q (u) given by [ID] 

1 

7j q 

with 



g q (u) = - u a ~ l e~ bu , (33) 



e* q = + (34) 

where [y] + = max(y, 0) and Z q is the normalization factor. Note that g q (u) reduces to the 
T distribution in the limit of q — > 1.0. Although the g~r distribution was useful for fs(u) 
of harmonic-oscillator systems subjected to finite bath [10], it does not work for fs{u) of 
double- well systems. This difference may be understood from a comparison between /s(w) 



for Ns = 1 of a double-well system shown in Fig. [I5[a) [or 16 a)] and its counterpart of a 
harmonic oscillator system shown in Fig. 9(a) of Ref. [TU]. Although the latter shows an 
exponential-like behavior with a monotonous decrease with increasing u, the former with 
a characteristic peak at u — 1.0 cannot be expressed by either the exponential, T, or q-T 
distribution. 



B. Analysis of stationary position distributions 

We have studied also the Ns dependence of stationary position distributions of p(Q) and 
P(Q), where Q denotes the position of a particle in the system and Q expresses the averaged 
position given by 

N S 

N s 



, N S 

q = f E^ ( 35 ) 



fe=l 



Figures 17|a) and 17 ^b) show p(Q) and -P(Q), respectively, obtained by DSs for various N s 
with N B = 100, T = 1.0 and c Q = 1.0. For N s = 1, we obtain p(Q) = P(Q) with the 
characteristic double-peaked structure. We note, however, that P{Q) is different from p{Q) 
for Ns > 1 for which P(Q) has a single-peaked structure despite the double-peaked p{Q). 
This is easily understood as follows: For example, in the case of Ns = 2, two particles in 
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FIG. 17: (Color online) Stationary distributions of (a) p(Q) as a function of particle position Q 
and (b) P(Q) as a function of the averaged position Q for various N$: N$ = 1 (dashed curve), 2 
(solid curve), 5 (dotted curve) and 10 (chain curve) with Nb = 100, T = 1.0 and c G = 1.0. Open 



circles in (b) express an analytical result obtained by Eq. (36) with N$ = 2 



the system mainly locate at Qk = 1.0 or Qk = —1.0 (k = 1,2) which yields the double- 
peaked distribution of fs(Q)- However, the averaged position of Q — (Qi + Q-i)/2 will be 
dominantly Q = 0.0, which leads to a single-peaked P(Q). The situation is the same also 
for N s > 2. 
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Theoretically P(Q) may be expressed by 



Ns / Ns \ 

P(Q) = [■ / 1] rfQfc exp [-/3V(Q k )} sIQ-N^qA. (36) 

k=l \ k=l / 

P(Q) numerically evaluated for Ns = 2 is plotted by open circles in Fig. 17 which are 
in good agreement with the solid curve expressing P(Q) obtained by DS. It is impossibly 



difficult to numerically evaluate Eq. (36) for N$ > 3. In the limit of Ns — > oo, P(Q) reduces 
to the Gaussian distribution according to the central-limit theorem. This trend is realized 
already in the case of Ns = 10 in Fig. [Tub). 



V. CONCLUDING REMARKS 



We have studied the properties of classical double-well systems coupled to finite bath, 
employing the (Ns + Nb) model [10] in which iVs-body system is coupled to iV^-body bath. 
Results obtained by DSs have shown the following: 

(i) Chaotic oscillations are induced in the double-well system coupled to finite bath in the 
absence of external forces for appropriate model parameters of c a , Nb, T, {u n } and Es , 

(ii) Among model parameters, fs(u) depends mainly on N s , c Q and T while fs(u) depends 
on Nb and T for N s < N B , 

(iii) fs(u) for Ns > 1 obtained by DSs may be phenomenological expressed by the T distri- 
bution, 

(iv) fs(u) for Ns = 1 with c Q ^ 0.0 cannot be described by either the exponential, T, or q-T 
distribution, although that with c = 0.0 follows the exponential distribution, and 

(v) The dissipation is not realized in the system energy for DSs at t = — 1000 with 
Ng = l- 100 and N B = 10 - 1000. 

The item (i) is in consistent with chaos in a closed classical double-well system driven by 
external forces [23J, although chaos is induced without external forces in our open classical 
double-well system. This is somewhat reminiscent of chaos induced by quantum noise in 
the absence of external force in closed quantum double- well systems [27] . Effects of induced 
chaos in the item (i) are not apparent in fs(u) because u (= us) is ensemble averaged over 10 
000 runs (realizations) with exponentially distributed initial system energies. Items (ii) and 
(v) are the same as in the harmonic-oscillator system coupled to finite bath [10]. The item 
(v) suggests that for the energy dissipation of system, we might need to adopt a much larger 
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Nb 1000) [26]. The item (iv) is in contrast to fs(u) for N$ = 1 in the open harmonic- 
oscillator system which may be approximately accounted for by the g~r distribution [TU] . It 
would be necessary and interesting to make a quantum extension of our study which is left 
as our future subject. 
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